package scu.maqiang.possion;

import java.util.Arrays;

import scu.maqiang.fes.FES3H81;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Possion3DH81 {

	public static void main(String[] args) {
        int N = 8, label = 10;
        double ep = 1.0 / N;
        
        //ʹ��lambda���ʽ�����ί��
        ScalarFunc regionFunc = (xy, llabel, param) -> {
            double xCell = xy[0] * N - (int)(xy[0] * N);
            double yCell = xy[1] * N - (int)(xy[1] * N);
            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
        };

		Mesh3H8 mesh3D = new Mesh3H8().cube(200, 200, 200);
        mesh3D.setDomainLabel(regionFunc, null, label);
        
        FES3H81 fs = new FES3H81(mesh3D);
        SRMatrix A = new SRMatrix(fs.GetNdof());

        ScalarFunc regionCoef = (xyz, relabel, param) -> relabel == label ? 1.0 : 1.0;
        ScalarFunc ConstCoef = (xyz, relabel, param) -> 1.0e0;
        ScalarFunc[] gcf = {regionCoef};
        gcf[0] = ConstCoef;
        
        ScalarFunc coefFunc = (xyz, relabel, param) -> 
                    2 * xyz[0] * (1 - xyz[0]) * xyz[1] * (1 - xyz[1])
                  + 2 * xyz[1] * (1 - xyz[1]) * xyz[2] * (1 - xyz[2])
                  + 2 * xyz[2] * (1 - xyz[2]) * xyz[0] * (1 - xyz[0]);
        fs.assembleStiff(gcf, null, BVPType.COMMON, A);
        //fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
        //System.out.println(A);
//        System.out.println(A.makeEntrySymmetric());
        double[] RHS = new double[fs.GetNdof()];
        //GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0;
        gcf[0] = ConstCoef;
        fs.assembleSource(gcf, null, BVPType.COMMON, RHS);
        A.makeEntrySymmetric();
        fs.applyBCZero(A, RHS, Direct.X, BCType.MBN, 1, 3, 5);
        
        //System.out.println(Arrays.toString(RHS));
        
//        System.out.println(A);

//        IterSSolver solver = new IterSSolver(1.0e-15, 1000000);
//        
        double[] x = new double[fs.GetNdof()];
//        CSRMatrix mat = new CSRMatrix(A);
//        solver.PCGSSOR(A, RHS, x, 1.5e0, 1);
//        Arrays.fill(x, 0.0);
//        solver.PCGSSOR_2(A, RHS, x, 1.5e0, 1);
//
//        
//        //System.out.println(A);
//        
//        //System.out.println(mat);
//        
//        Arrays.fill(x, 0.0);
//        solver.CG(mat, RHS, x, 1);
//        
//        Arrays.fill(x, 0.0);
//        solver.CG(A, RHS, x, 1);
        
        Arrays.fill(x, 0.0);
        NewIterSSolver newSolver = new NewIterSSolver(A, 1.0e-15, 1000000);
        newSolver.PCGSSOR(RHS, x, 1.5e0, 1);

        Arrays.fill(x, 0.0);
        //newSolver.GMRES(RHS, x, 200);
        newSolver.GMRES_SSOR(RHS, x, 1.5, 30, 1);
        mesh3D.toTecplot("Solutions3DH8.dat", x);
        
        double[][] gradX = new double[3][fs.GetNdof()];
        fs.computeGradient2(x, gradX);
        mesh3D.toTecplot("Solutions.dat", x);
        mesh3D.toTecplot("SolutionsGradient.dat", gradX);
        
        System.out.println("Gradient: ");
        System.out.println("Dx min: " + MVO.min(gradX[0]) + "\tmax: " + MVO.max(gradX[0]));
        System.out.println("Dy min: " + MVO.min(gradX[1]) + "\tmax: " + MVO.max(gradX[1]));
        System.out.println("Dz min: " + MVO.min(gradX[2]) + "\tmax: " + MVO.max(gradX[2]));
        
//        ArrayList<Integer>  li = new ArrayList<Integer>();
//        li.ensureCapacity(100);
//        for(int i = 0; i < 100; i++) {
//        	li.add(null);
//        }
//        li.add(1, 2);
//        System.out.println(li.toString());
        
//        Mesh2DT3 mesh6 = mg.square2D(2, 2, T3Type.Left);
//        mesh6.displayBoundarys();
//        mesh6.extractBoundarysFromElements();
//        mesh6.displayBoundarys();
	}
	
	/**
	 * Possion����Dirichlet���������
	 * <p>
	 * &part;<sub>i</sub>(a<sub>ij</sub>&part;<sub>j</sub>u) = f &#915;
	 * 
	 * @param kxxFunc
	 * @param paramKxx
	 * @param sourceFunc
	 * @param paramSource
	 * @param tp
	 * @param DirichletBoundaryLabel
	 * @param fixedFunc
	 * @param paramFixed
	 * @param reportFile
	 */
	public void DirichletProblem(BVPType tp, ScalarFunc kxxFunc, double[] paramKxx, ScalarFunc sourceFunc,
			double[] paramSource, int[] DirichletBoundaryLabel, ScalarFunc fixedFunc, double[] paramFixed,
			String reportFile) {
		fs.assembleStiff(new ScalarFunc[] { kxxFunc }, paramKxx, tp, A);
		if (sourceFunc != null) {
			fs.assembleSource(new ScalarFunc[] { sourceFunc }, paramSource, tp, RHS);
		}
		
		fs.applyBC_MBN(A, Direct.All, DirichletBoundaryLabel);

		if (fixedFunc == null) {
			fs.applyBC_MBN(RHS, Direct.All, 0.0, DirichletBoundaryLabel);
		} else {
			fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[] { fixedFunc }, paramFixed, DirichletBoundaryLabel);
		}
		System.out.println("OK!");
		solveAndComputeGradientAndExport(A, RHS, x, reportFile);
	}
	
	private void solveAndComputeGradientAndExport(SRMatrix A, double[] RHS, double[] x, String reportFile) {
		NewIterSSolver solver = new NewIterSSolver(A, true);
		solver.PCGSSSOR(RHS, x, 1.50, 2);
		fs.computeGradient(x, gradient);
		System.out.println("Solution: ");
		System.out.println("x0 min:" + MVO.min(x) + "\t" + "max: " + MVO.max(x));
		System.out.println("Solution Gradient: ");
		System.out.println("Dx min:" + MVO.min(gradient[0]) + "\t" + "max: " + MVO.max(gradient[0]));
		System.out.println("Dy min:" + MVO.min(gradient[1]) + "\t" + "max: " + MVO.max(gradient[1]));
		System.out.println("Dz min:" + MVO.min(gradient[2]) + "\t" + "max: " + MVO.max(gradient[2]));
		double[][] mat = new double[][] { x, gradient[0], gradient[1], gradient[2]};
		if (reportFile != null) {
			mesh.toTecplot(reportFile, mat);
		}
	}
	
	
	public double[] exportSolution() {
		return x.clone();
	}
	Mesh3H8 mesh;
	FES3H81 fs;
	SRMatrix A;
	double[] RHS;
	double[] x;
	double[][] gradient;

}
